clear all
set more off
set mem 10000000
set matsize 10000
version 15

**************************************************
*** RDROBUST first stage, nighttime brightness ***
**************************************************

** Set file paths
do "$path_code/paths.do"

** Set graph scheme
cd "$path/code/analyze"
set scheme fb, perm

****************************************************************** 
****************************************************************** 

** 1. Main lights RD: RDROBUST, tri kernel, default SE, STFEs, pre lights ctrls 
**                    (1998-2004), drop lights_diff>=20, drop pop_mismatch20
{
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in RD sample
gen in_fs_sample = vplan4<11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==1 & pop_non_zero==1 
keep if in_fs_sample==1
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==32 // Kerala, only 3 villages
tab st_code, gen(STFE)	
drop STFE1 // to avoid collinearity

	// Create district FEs (since RD robust doesn't let you pass them through)
gen stdtFE = stdt
tab stdtFE state     if inlist(stdtFE,62,63,64,70,91,180,181,202,308,309,320,490,493,494,495,504,505,508)	
replace stdtFE = 999 if inlist(stdtFE,62,63,64,70,91,180,181,202,308,309,320,490,493,494,495,504,505,508)	
	// one catch-all district FE for districts with so few in-sample villages that they break rdrobust
tab stdtFE, gen(DTFE)	
drop DTFE1 // to avoid collinearity

	// Create block groups (for clustering)
egen stdtbk = group(stdt bk_code)

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)


  // Create variables to store regresson results
gen fs_step = .
gen yvar = ""
gen ifs = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen ndist = . 
gen ymean = .
gen ftag = ""
local folder = "RDROBUST plots fs lights"
local h_wide = 200


	// 1.1 Lights brightness regressions, preferred
local fs_step = 1
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1
local row = 0

foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

	// Define outcome variable
	local yvar = "lights_max`y'_hat"
	local title = "Preferred specification: `y' Brightness"
	local ftag = "preferred"

	// Remove controls used to project outcome variable
	local controls = "`CONTROLS'"
	if regexm("`yvar'","2007")==1 {
		local controls = subinstr("`controls'"," lights_max2005","",1)
	}
	else if regexm("`yvar'","2006")==1 {
		local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2005")==1 {
		local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2004")==1 {
		local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2003")==1 {
		local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2002")==1 {
		local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	
	// Run first-stage regresssion
	rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

	// Generate in-sample indicator and store bandwidths
	cap drop temp_in_reg
	qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
	qui count if temp_in_reg
	assert r(N)==(e(N_h_l)+e(N_h_r))
	local h_l = e(h_l)
	local h_r = e(h_r)

	// Store results
	local row = `row' + 1
	qui replace fs_step = `fs_step' in `row'
	qui replace yvar = "`yvar'" in `row'
	qui replace ifs = "`ifs'" in `row'
	qui replace control = "`controls'" in `row'
	qui replace fe = "`fe'" in `row'
	qui replace kernel = e(kernel) in `row'
	qui replace bwmethod = e(bwselect) in `row'
	qui replace vce = "`vce'" in `row'
	qui replace polynomial_order = e(p) in `row'
	qui replace beta_conv = e(tau_cl) in `row'
	qui replace beta_robust = e(tau_bc) in `row'
	qui replace se_conv = e(se_tau_cl) in `row'
	qui replace se_robust = e(se_tau_rb) in `row'
	qui replace pval_conv = e(pv_cl) in `row'
	qui replace pval_robust = e(pv_rb) in `row'
	qui replace lci_conv = e(ci_l_cl) in `row'
	qui replace uci_conv = e(ci_r_cl) in `row'
	qui replace lci_robust = e(ci_l_rb) in `row'
	qui replace uci_robust = e(ci_r_rb) in `row'
	qui replace bw_lo = e(h_l) in `row'
	qui replace bw_hi = e(h_r) in `row'
	qui replace nobs_orig = e(N) in `row'
	qui replace nobs_left = e(N_h_l) in `row'
	qui replace nobs_right = e(N_h_r) in `row'
	qui unique stdt if temp_in_reg==1
	qui replace ndist = r(unique) in `row'
	qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
	qui replace ymean = r(mean) in `row'
	qui replace ftag = "`ftag'" in `row'

	// Residualize `yvar' for regression sample
	qui reg `yvar' `controls' `fe' if temp_in_reg==1
	cap drop y_resid_sample
	predict y_resid_sample, residuals
	
	// Residualize `yvar' for [100,500] bandwidth
	qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
	cap drop y_resid_full
	predict y_resid_full, residuals
	
	// RD plots for regression sample (variable bandwidth) 
	foreach bins in 10 20 40 {
		if `h_l'<=75 {
			local xlab = "250 275 300 325 350"
		}
		else if inrange(`h_l',76,125) {
			local xlab = "200 250 300 350 400"
		}
		else if inrange(`h_l',126,175) {
			local xlab = "150 200 250 300 350 400 450"
		}
		else {
			local xlab = ""
		}
		rdplot y_resid_sample tot_p if temp_in_reg==1, ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`y' brightness residuals", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(`xlab', labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
	}

	// RD plots for regression sample (constant bandwidth) 
	foreach bins in 10 20 40 {
		if `bins'==40 {
			local ymarks = "-.4 -.2 0 .2 .4"
		}
		else if `bins'==20 {
			local ymarks = "-.3 -.2 -.1 0 .1 .2 .3"
		}
		else if `bins'==10 {
			local ymarks = "-.2 -.1 0 .1 .2"
		}
		rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`y' brightness residuals", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(`ymarks',nogrid angle(0) labsize(medlarge)) ///
			xlabel(, labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
		local ymarks = ""
	}			
}



	// 1.2 Sensitivity to number of projections years
local fs_step = 2
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1

foreach y in 2011 {

	foreach pyrs in 0 1 {

		// Define outcome variable
		if `pyrs'==0 {
			local yvar = "lights_max`y'"
			local title = "Raw `y' Brightness"
			local ftag = "raw"
		}
		if `pyrs'==1 {
			local yvar = "lights_max`y'_1_hat"
			local title = "Projected `y' Brightness, Adjacent Years"
			local ftag = "adjyrs"
		}

		// Remove controls used to project outcome variable
		local controls = "`CONTROLS'"
		if regexm("`yvar'","2007")==1 {
			local controls = subinstr("`controls'"," lights_max2005","",1)
		}
		else if regexm("`yvar'","2006")==1 {
			local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2005")==1 {
			local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2004")==1 {
			local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2003")==1 {
			local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2002")==1 {
			local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
		qui count if temp_in_reg
		assert r(N)==(e(N_h_l)+e(N_h_r))
		local h_l = e(h_l)
		local h_r = e(h_r)
		
		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'

		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// Residualize `yvar' for [100,500] bandwidth
		qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
		cap drop y_resid_full
		predict y_resid_full, residuals
		
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
		}

		// RD plots for regression sample (constant bandwidth) 
		foreach bins in 10 20 40 {
			rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(, labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
		}
	}		
}



	// 1.3 Sensitivity to different lights variables
local fs_step = 3
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1

foreach y in 2011 {

	foreach reg in 1 2 3 4 5 {

		// Define outcome variable
		local pre1  = substr("0" + string(`y'-2000-1),-2,2)
		local pre2  = substr("0" + string(`y'-2000-2),-2,2)
		local post1 = substr("0" + string(`y'-2000+1),-2,2)
		local post2 = substr("0" + string(`y'-2000+2),-2,2)
		if `reg'==1 {
			local yvar = "lights_mean`y'_hat"
			local title = "Mean `y' Brightness, Projected"
			local ftag = "meanhat"
		}
		if `reg'==2 {
			local yvar = "post_avg_max_`pre1'_`post1'"
			local title = "Maximum `y' Brightness, 3-Year Average"
			local ftag = "max3yravg"
		}
		if `reg'==3 {
			local yvar = "post_avg_mean_`pre1'_`post1'"
			local title = "Mean `y' Brightness, 3-Year Average"
			local ftag = "mean3yravg"
		}
		if `reg'==4 {
			local yvar = "post_avg_max_`pre2'_`post2'"
			local title = "Maximum `y' Brightness, 5-Year Average"
			local ftag = "max5yravg"
		}
		if `reg'==5 {
			local yvar = "post_avg_mean_`pre2'_`post2'"
			local title = "Mean `y' Brightness, 5-Year Average"
			local ftag = "mean5yravg"
		}

		// Remove controls used to project outcome variable
		local controls = "`CONTROLS'"
		if regexm("`yvar'","2007")==1 {
			local controls = subinstr("`controls'"," lights_max2005","",1)
		}
		else if regexm("`yvar'","2006")==1 {
			local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2005")==1 {
			local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2004")==1 {
			local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2003")==1 {
			local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2002")==1 {
			local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
		qui count if temp_in_reg
		assert r(N)==(e(N_h_l)+e(N_h_r))
		local h_l = e(h_l)
		local h_r = e(h_r)
		
		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'

		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// Residualize `yvar' for [100,500] bandwidth
		qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
		cap drop y_resid_full
		predict y_resid_full, residuals
		
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
		}

		// RD plots for regression sample (constant bandwidth) 
		foreach bins in 10 20 40 {
			rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(, labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
		}
	}		
}



	// 1.4 Sensitivity to average vs. stable lights
local fs_step = 4
local ifs = "pop_mismatch20==0 & lights_diff<20"
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1

foreach y in 2011 {

	foreach reg in 1 2  {

		// Define outcome variable
		if `reg'==1 {
			local yvar = "lights_sta_max`y'_hat"
			local CONTROLS = "lights_sta_max1998 lights_sta_max1999 lights_sta_max2000 lights_sta_max2001 lights_sta_max2002 lights_sta_max2003 lights_sta_max2004 lights_sta_max2005" 
			local title = "`y' Brightness, Stable Lights"
			local ftag = "stablelights"
		}
		if `reg'==2 {
			local yvar = "lights_pct_max`y'_hat"
			local CONTROLS = "lights_pct_max1998 lights_pct_max1999 lights_pct_max2000 lights_pct_max2001 lights_pct_max2002 lights_pct_max2003 lights_pct_max2004 lights_pct_max2005" 
			local title = "`y' Brightness, Percent Lights"
			local ftag = "percentlights"
		}
		
		// Remove controls used to project outcome variable
		local controls = "`CONTROLS'"
		if regexm("`yvar'","2007")==1 {
			local controls = subinstr("`controls'"," lights_max2005","",1)
		}
		else if regexm("`yvar'","2006")==1 {
			local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2005")==1 {
			local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2004")==1 {
			local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2003")==1 {
			local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2002")==1 {
			local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
		qui count if temp_in_reg
		assert r(N)==(e(N_h_l)+e(N_h_r))
		local h_l = e(h_l)
		local h_r = e(h_r)
		
		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'

		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// Residualize `yvar' for [100,500] bandwidth
		qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
		cap drop y_resid_full
		predict y_resid_full, residuals
		
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
		}

		// RD plots for regression sample (constant bandwidth) 
		foreach bins in 10 20 40 {
			rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(, labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
		}
	}		
}



	// 1.5 Sensitivity to including population mismatches
local fs_step = 5
local ifs = "lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1

foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

	// Define outcome variable
	local yvar = "lights_max`y'_hat"
	local title = "`y' Brightness, Including Population Mismatches"
	local ftag = "inclpop20"
	
	// Remove controls used to project outcome variable
	local controls = "`CONTROLS'"
	if regexm("`yvar'","2007")==1 {
		local controls = subinstr("`controls'"," lights_max2005","",1)
	}
	else if regexm("`yvar'","2006")==1 {
		local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2005")==1 {
		local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2004")==1 {
		local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2003")==1 {
		local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2002")==1 {
		local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	
	// Run first-stage regresssion
	rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

	// Generate in-sample indicator and store bandwidths
	cap drop temp_in_reg
	qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
	qui count if temp_in_reg
	assert r(N)==(e(N_h_l)+e(N_h_r))
	local h_l = e(h_l)
	local h_r = e(h_r)
	
	// Store results
	local row = `row' + 1
	qui replace fs_step = `fs_step' in `row'
	qui replace yvar = "`yvar'" in `row'
	qui replace ifs = "`ifs'" in `row'
	qui replace control = "`controls'" in `row'
	qui replace fe = "`fe'" in `row'
	qui replace kernel = e(kernel) in `row'
	qui replace bwmethod = e(bwselect) in `row'
	qui replace vce = "`vce'" in `row'
	qui replace polynomial_order = e(p) in `row'
	qui replace beta_conv = e(tau_cl) in `row'
	qui replace beta_robust = e(tau_bc) in `row'
	qui replace se_conv = e(se_tau_cl) in `row'
	qui replace se_robust = e(se_tau_rb) in `row'
	qui replace pval_conv = e(pv_cl) in `row'
	qui replace pval_robust = e(pv_rb) in `row'
	qui replace lci_conv = e(ci_l_cl) in `row'
	qui replace uci_conv = e(ci_r_cl) in `row'
	qui replace lci_robust = e(ci_l_rb) in `row'
	qui replace uci_robust = e(ci_r_rb) in `row'
	qui replace bw_lo = e(h_l) in `row'
	qui replace bw_hi = e(h_r) in `row'
	qui replace nobs_orig = e(N) in `row'
	qui replace nobs_left = e(N_h_l) in `row'
	qui replace nobs_right = e(N_h_r) in `row'
	qui unique stdt if temp_in_reg==1
	qui replace ndist = r(unique) in `row'
	qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
	qui replace ymean = r(mean) in `row'
	qui replace ftag = "`ftag'" in `row'

	// Residualize `yvar' for regression sample
	qui reg `yvar' `controls' `fe' if temp_in_reg==1
	cap drop y_resid_sample
	predict y_resid_sample, residuals
	
	// Residualize `yvar' for [100,500] bandwidth
	qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
	cap drop y_resid_full
	predict y_resid_full, residuals
	
	// RD plots for regression sample (variable bandwidth) 
	foreach bins in 10 20 40 {
		if `h_l'<=75 {
			local xlab = "250 275 300 325 350"
		}
		else if inrange(`h_l',76,125) {
			local xlab = "200 250 300 350 400"
		}
		else if inrange(`h_l',126,175) {
			local xlab = "150 200 250 300 350 400 450"
		}
		else {
			local xlab = ""
		}
		rdplot y_resid_sample tot_p if temp_in_reg==1, ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`y' brightness residuals", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(`xlab', labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
	}

	// RD plots for regression sample (constant bandwidth) 
	foreach bins in 10 20 40 {
		rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`y' brightness residuals", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(, labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
	}		
}



	// 1.6 Sensitivity to lights_diff outliers
local fs_step = 6
local ifs_base = "pop_mismatch20==0"
local CONTOOLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1

foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

	foreach ld in 10 15 25 30 35 40 45 50  {

		// Define outcome variable
		local yvar = "lights_max`y'_hat"
		
		// Define threshould for dropping lights_diff outliers
		local ifs = "`ifs_base' & lights_diff<`ld'"
		local title = "`y' Brightness, Dropping >`ld' Brightness Changes"
		local ftag = "ldiff`ld'"
		
		// Remove controls used to project outcome variable
		local controls = "`CONTROLS'"
		if regexm("`yvar'","2007")==1 {
			local controls = subinstr("`controls'"," lights_max2005","",1)
		}
		else if regexm("`yvar'","2006")==1 {
			local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2005")==1 {
			local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2004")==1 {
			local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2003")==1 {
			local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2002")==1 {
			local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
		qui count if temp_in_reg
		assert r(N)==(e(N_h_l)+e(N_h_r))
		local h_l = e(h_l)
		local h_r = e(h_r)
		
		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'

		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// Residualize `yvar' for [100,500] bandwidth
		qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
		cap drop y_resid_full
		predict y_resid_full, residuals
		
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
		}

		// RD plots for regression sample (constant bandwidth) 
		foreach bins in 10 20 40 {
			rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(, labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
		}		
	}
}



	// 1.7 Sensitivity to different bandwidth calculations
local fs_step = 7
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local vce = ""
local poly = 1

foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

	foreach reg in 1 2 3 4 5 {

		// Define outcome variable
		local yvar = "lights_max`y'_hat"

		// Define bandwidth method
		if `reg'==1 {
			local bwmethod = "msetwo"
			local title = "`y' Brightness, Bandwidth Method `bwmethod'"
			local ftag = "BW`bwmethod'"
		}
		if `reg'==2 {
			local bwmethod = "msesum"
			local title = "`y' Brightness, Bandwidth Method `bwmethod'"
			local ftag = "BW`bwmethod'"
		}
		if `reg'==3 {
			local bwmethod = "cerrd"
			local title = "`y' Brightness, Bandwidth Method `bwmethod'"
			local ftag = "BW`bwmethod'"
		}
		if `reg'==4 {
			local bwmethod = "certwo"
			local title = "`y' Brightness, Bandwidth Method `bwmethod'"
			local ftag = "BW`bwmethod'"
		}
		if `reg'==5 {
			local bwmethod = "cersum"
			local title = "`y' Brightness, Bandwidth Method `bwmethod'"
			local ftag = "BW`bwmethod'"
		}

		
		// Remove controls used to project outcome variable
		local controls = "`CONTROLS'"
		if regexm("`yvar'","2007")==1 {
			local controls = subinstr("`controls'"," lights_max2005","",1)
		}
		else if regexm("`yvar'","2006")==1 {
			local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2005")==1 {
			local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2004")==1 {
			local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2003")==1 {
			local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2002")==1 {
			local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
		qui count if temp_in_reg
		assert r(N)==(e(N_h_l)+e(N_h_r))
		local h_l = e(h_l)
		local h_r = e(h_r)
		
		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'

		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
		}
	}		
}



	// 1.8 Sensitivity to 2nd-order polynomial and differet kernels
local fs_step = 8
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local bwmethod = "mserd"
local vce = ""

foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

	foreach reg in 1 2 3 4 5 {

		// Define outcome variable
		local yvar = "lights_max`y'_hat"

		// Define kernel and polynomial combo
		if `reg'==1 {
			local kernel = "tri"
			local poly = 2
			local title = "`y' Brightness, Triangular Kernel, Quadratic"
			local ftag = "K`kernel'P`poly'"
		}
		if `reg'==2 {
			local kernel = "epa"
			local poly = 1
			local title = "`y' Brightness, Epanechnikov Kernel"
			local ftag = "K`kernel'P`poly'"
		}
		if `reg'==3 {
			local kernel = "epa"
			local poly = 2
			local title = "`y' Brightness, Epanechnikov Kernel, Quadratic"
			local ftag = "K`kernel'P`poly'"
		}
		if `reg'==4 {
			local kernel = "uni"
			local poly = 1
			local title = "`y' Brightness, Uniform Kernel"
			local ftag = "K`kernel'P`poly'"
		}
		if `reg'==5 {
			local kernel = "uni"
			local poly = 2
			local title = "`y' Brightness, Uniform Kernel, Quadratic"
			local ftag = "K`kernel'P`poly'"
		}

		
		// Remove controls used to project outcome variable
		local controls = "`CONTROLS'"
		if regexm("`yvar'","2007")==1 {
			local controls = subinstr("`controls'"," lights_max2005","",1)
		}
		else if regexm("`yvar'","2006")==1 {
			local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2005")==1 {
			local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2004")==1 {
			local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2003")==1 {
			local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2002")==1 {
			local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
		qui count if temp_in_reg
		assert r(N)==(e(N_h_l)+e(N_h_r))
		local h_l = e(h_l)
		local h_r = e(h_r)
		
		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'

		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// Residualize `yvar' for [100,500] bandwidth
		qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
		cap drop y_resid_full
		predict y_resid_full, residuals
		
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
		}

		// RD plots for regression sample (constant bandwidth) 
		foreach bins in 10 20 40 {
			rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(, labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
		}
	}		
}



	// 1.9 Sensitivity to fixed effects
local fs_step = 9
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1

foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

	foreach reg in 1 2 {

		// Define outcome variable
		local yvar = "lights_max`y'_hat"

		// Define fixed effects
		if `reg'==1 {
			local fe = ""
			local title = "`y' Brightness, No Fixed Effects"
			local ftag = "nofes"
		}
		if `reg'==2 {
			local fe = "DTFE*"
			local title = "`y' Brightness, District Fixed Effects"
			local ftag = "dtfes"
		}

		// Remove controls used to project outcome variable
		local controls = "`CONTROLS'"
		if regexm("`yvar'","2007")==1 {
			local controls = subinstr("`controls'"," lights_max2005","",1)
		}
		else if regexm("`yvar'","2006")==1 {
			local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2005")==1 {
			local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2004")==1 {
			local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2003")==1 {
			local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2002")==1 {
			local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
		qui count if temp_in_reg
		assert r(N)==(e(N_h_l)+e(N_h_r))
		local h_l = e(h_l)
		local h_r = e(h_r)
		
		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'

		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// Residualize `yvar' for [100,500] bandwidth
		qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
		cap drop y_resid_full
		predict y_resid_full, residuals
		
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
		}

		// RD plots for regression sample (constant bandwidth) 
		foreach bins in 10 20 40 {
			rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(, labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
		}
	}		
}



	// 1.10 Sensitivity to pre controls
local fs_step = 10
local ifs = "pop_mismatch20==0 & lights_diff<20"
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1

foreach y in 2005 2006 2007 2008 2009 2010 2011 2012 {

	foreach reg in 1 2 3 4 5 6 7 8 9 10 {

		// Define outcome variable
		local yvar = "lights_max`y'_hat"

		// Define fixed effects
		if `reg'==1 {
			local CONTROLS = "lights_max2001" 
			local title = "`y' Brightness, 2001 Lights Control"
			local ftag = "ctrl01"
		}
		if `reg'==2 {
			local CONTROLS = "lights_max2000 lights_max2001 lights_max2002" 
			local title = "`y' Brightness, 2000-2002 Lights Controls"
			local ftag = "ctrl0002"
		}
		if `reg'==3 {
			local CONTROLS = "lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003" 
			local title = "`y' Brightness, 1999-2003 Lights Controls"
			local ftag = "ctrl9903"
		}
		if `reg'==4 {
			local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004" 
			local title = "`y' Brightness, 1998-2004 Lights Controls"
			local ftag = "ctrl9804"
		}
		if `reg'==5 {
			local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 vdpr_tra_d_app_pr_01" 
			local title = "`y' Brightness, Control: 2001 presence of road"
			local ftag = "ctrlroad"
		}
		if `reg'==6 {
			local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 pct_irr_01" 
			local title = "`y' Brightness, Control: 2001 share of land irrigated"
			local ftag = "ctrlpctirr"
		}
		if `reg'==7 {
			local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 lit_p_01" 
			local title = "`y' Brightness, Control: 2001 literacy"
			local ftag = "ctrllit"
		}
		if `reg'==8 {
			local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 pct_scst_01" 
			local title = "`y' Brightness, Control: 2001 share SC/ST"
			local ftag = "ctrlscst"
		}
		if `reg'==9 {
			local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 vd_pwr_d_any_01" 
			local title = "`y' Brightness, Control: 2001 power in any sector"
			local ftag = "ctrlpwrany"
		}
		if `reg'==10 {
			local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005 vd_pwr_d_all_01" 
			local title = "`y' Brightness, Control: 2001 power in 3 sectors"
			local ftag = "ctrlpwrall"
		}

		// Remove controls used to project outcome variable
		local controls = "`CONTROLS'"
		if regexm("`yvar'","2007")==1 {
			local controls = subinstr("`controls'"," lights_max2005","",1)
		}
		else if regexm("`yvar'","2006")==1 {
			local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2005")==1 {
			local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2004")==1 {
			local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2003")==1 {
			local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2002")==1 {
			local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
		*qui count if temp_in_reg
		*assert r(N)==(e(N_h_l)+e(N_h_r)) // added controls are occasionaly missing
		local h_l = e(h_l)
		local h_r = e(h_r)
		
		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'

		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// Residualize `yvar' for [100,500] bandwidth
		qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
		cap drop y_resid_full
		predict y_resid_full, residuals
		
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
		}

		// RD plots for regression sample (constant bandwidth) 
		foreach bins in 10 20 40 {
			rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(, labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
		}
	}		
}



	// 1.11 Sensitivity to alternative standard errors
local fs_step = 12
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local poly = 1

foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

	foreach reg in 1 2 3 4 5 6 7 {

		// Define outcome variable
		local yvar = "lights_max`y'_hat"

		// Define vce method
		if `reg'==1 {
			local vce = "vce(nn 5)"
			local title = "`y' Brightness, Nearest-Neighbor SEs (min 5)"
			local ftag = "VCEnn5"
		}
		if `reg'==2 {
			local vce = "vce(nncluster stdt)"
			local title = "`y' Brightness, NN-Clustered SEs (District)"
			local ftag = "VCEnnclDT"
		}
		if `reg'==3 {
			local vce = "vce(nncluster stdtbk)"
			local title = "`y' Brightness, NN-Clustered SEs (Block)"
			local ftag = "VCEnnclBK"
		}
		if `reg'==4 {
			local vce = "vce(nncluster stdt 5)"
			local title = "`y' Brightness, NN-Clustered SEs (District, min 5)"
			local ftag = "VCEnnclDT5"
		}
		if `reg'==5 {
			local vce = "vce(nncluster stdtbk 5)"
			local title = "`y' Brightness, NN-Clustered SEs (Block, min 5)"
			local ftag = "VCEnnclBK5"
		}
		if `reg'==6 {
			local vce = "vce(cluster stdt)"
			local title = "`y' Brightness, Clustered SEs (District)"
			local ftag = "VCEclDT"
		}
		if `reg'==7 {
			local vce = "vce(cluster stdtbk)"
			local title = "`y' Brightness, Clustered SEs (Block)"
			local ftag = "VCEclBK"
		}
		
		// Remove controls used to project outcome variable
		local controls = "`CONTROLS'"
		if regexm("`yvar'","2007")==1 {
			local controls = subinstr("`controls'"," lights_max2005","",1)
		}
		else if regexm("`yvar'","2006")==1 {
			local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2005")==1 {
			local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2004")==1 {
			local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2003")==1 {
			local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		else if regexm("`yvar'","2002")==1 {
			local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
		}
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
		qui count if temp_in_reg
		assert r(N)==(e(N_h_l)+e(N_h_r))
		local h_l = e(h_l)
		local h_r = e(h_r)
		
		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = e(bwselect) in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'

	}		
}


	// Save results
keep fs_step-ftag
dropmiss, obs force
compress
save "$results/RDROBUST_fs_lights.dta", replace


}

****************************************************************** 
****************************************************************** 

** 2. Falsification tests for main lights RD
{

	// 2.1 11th Plan, single habitation
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in FALSIFICATION sample: 11th Plan, single habitation
gen in_fs_sampleFALS = vplan4==11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==1 & pop_non_zero==1 
keep if in_fs_sampleFALS==1
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==14 // Manipur has only 7 between 100 and 500 people
drop if st_code==32 // Kerala, only 3 villages
drop if st_code==33 // Tamil Nadu has lots of village, but very few beween 200 and 400 people
tab st_code, gen(STFE)	
drop STFE1 // to avoid collinearity

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)

  // Create variables to store regresson results
gen fs_step = .
gen yvar = ""
gen ifs = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen ndist = . 
gen ymean = .
gen ftag = ""
gen fals = ""
local folder = "RDROBUST plots fs lights"
local h_wide = 200

local fs_step = 13
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1
local fals = "11th Plan, single hab"
local row = 0

foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

	// Define outcome variable
	local yvar = "lights_max`y'_hat"
	local title = "`y' Brightness, 11th Plan, Single Habitation"
	local ftag = "FALS_11th_singH"

	// Remove controls used to project outcome variable
	local controls = "`CONTROLS'"
	if regexm("`yvar'","2007")==1 {
		local controls = subinstr("`controls'"," lights_max2005","",1)
	}
	else if regexm("`yvar'","2006")==1 {
		local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2005")==1 {
		local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2004")==1 {
		local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2003")==1 {
		local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2002")==1 {
		local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}

	// Run first-stage regresssion
	rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

	// Generate in-sample indicator and store bandwidths
	cap drop temp_in_reg
	qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
	qui count if temp_in_reg
	assert r(N)==(e(N_h_l)+e(N_h_r))
	local h_l = e(h_l)
	local h_r = e(h_r)

	// Store results
	local row = `row' + 1
	qui replace fs_step = `fs_step' in `row'
	qui replace yvar = "`yvar'" in `row'
	qui replace ifs = "`ifs'" in `row'
	qui replace control = "`controls'" in `row'
	qui replace fe = "`fe'" in `row'
	qui replace kernel = e(kernel) in `row'
	qui replace bwmethod = e(bwselect) in `row'
	qui replace vce = "`vce'" in `row'
	qui replace polynomial_order = e(p) in `row'
	qui replace beta_conv = e(tau_cl) in `row'
	qui replace beta_robust = e(tau_bc) in `row'
	qui replace se_conv = e(se_tau_cl) in `row'
	qui replace se_robust = e(se_tau_rb) in `row'
	qui replace pval_conv = e(pv_cl) in `row'
	qui replace pval_robust = e(pv_rb) in `row'
	qui replace lci_conv = e(ci_l_cl) in `row'
	qui replace uci_conv = e(ci_r_cl) in `row'
	qui replace lci_robust = e(ci_l_rb) in `row'
	qui replace uci_robust = e(ci_r_rb) in `row'
	qui replace bw_lo = e(h_l) in `row'
	qui replace bw_hi = e(h_r) in `row'
	qui replace nobs_orig = e(N) in `row'
	qui replace nobs_left = e(N_h_l) in `row'
	qui replace nobs_right = e(N_h_r) in `row'
	qui unique stdt if temp_in_reg==1
	qui replace ndist = r(unique) in `row'
	qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
	qui replace ymean = r(mean) in `row'
	qui replace ftag = "`ftag'" in `row'
	qui replace fals = "`fals'" in `row'

	// Residualize `yvar' for regression sample
	qui reg `yvar' `controls' `fe' if temp_in_reg==1
	cap drop y_resid_sample
	predict y_resid_sample, residuals
	
	// Residualize `yvar' for [100,500] bandwidth
	qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
	cap drop y_resid_full
	predict y_resid_full, residuals
	
	// RD plots for regression sample (variable bandwidth) 
	foreach bins in 10 20 40 {
		if `h_l'<=75 {
			local xlab = "250 275 300 325 350"
		}
		else if inrange(`h_l',76,125) {
			local xlab = "200 250 300 350 400"
		}
		else if inrange(`h_l',126,175) {
			local xlab = "150 200 250 300 350 400 450"
		}
		else {
			local xlab = ""
		}
		rdplot y_resid_sample tot_p if temp_in_reg==1, ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`y' brightness residuals", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(`xlab', labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
	}

	// RD plots for regression sample (constant bandwidth) 
	foreach bins in 10 20 40 {
		rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`y' brightness residuals", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(, labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
	}			
}

	// Store results in temp file
keep fs_step-ftag fals
dropmiss, obs force
tempfile fals1
save `fals1'



	// 2.2 10th Plan, multi habitation
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in FALSIFICATION sample: 10th Plan, multi habitation
gen in_fs_sampleFALS = vplan4<11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==0 & pop_non_zero==1 
keep if in_fs_sampleFALS==1
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==6 // Haryana has only 1 village in bandwidth
drop if st_code==32 // Kerala has no villages in bandwidth
tab st_code, gen(STFE)	
drop STFE1 // to avoid collinearity

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)

  // Create variables to store regresson results
gen fs_step = .
gen yvar = ""
gen ifs = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen ndist = . 
gen ymean = .
gen ftag = ""
gen fals = ""
local folder = "RDROBUST plots fs lights"
local h_wide = 200

local fs_step = 13
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1
local fals = "10th Plan, multi hab"
local row = 0

foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

	// Define outcome variable
	local yvar = "lights_max`y'_hat"
	local title = "`y' Brightness, 10th Plan, Multiple Habitations"
	local ftag = "FALS_10th_multiH"

	// Remove controls used to project outcome variable
	local controls = "`CONTROLS'"
	if regexm("`yvar'","2007")==1 {
		local controls = subinstr("`controls'"," lights_max2005","",1)
	}
	else if regexm("`yvar'","2006")==1 {
		local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2005")==1 {
		local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2004")==1 {
		local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2003")==1 {
		local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2002")==1 {
		local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}

	// Run first-stage regresssion
	rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

	// Generate in-sample indicator and store bandwidths
	cap drop temp_in_reg
	qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
	qui count if temp_in_reg
	assert r(N)==(e(N_h_l)+e(N_h_r))
	local h_l = e(h_l)
	local h_r = e(h_r)

	// Store results
	local row = `row' + 1
	qui replace fs_step = `fs_step' in `row'
	qui replace yvar = "`yvar'" in `row'
	qui replace ifs = "`ifs'" in `row'
	qui replace control = "`controls'" in `row'
	qui replace fe = "`fe'" in `row'
	qui replace kernel = e(kernel) in `row'
	qui replace bwmethod = e(bwselect) in `row'
	qui replace vce = "`vce'" in `row'
	qui replace polynomial_order = e(p) in `row'
	qui replace beta_conv = e(tau_cl) in `row'
	qui replace beta_robust = e(tau_bc) in `row'
	qui replace se_conv = e(se_tau_cl) in `row'
	qui replace se_robust = e(se_tau_rb) in `row'
	qui replace pval_conv = e(pv_cl) in `row'
	qui replace pval_robust = e(pv_rb) in `row'
	qui replace lci_conv = e(ci_l_cl) in `row'
	qui replace uci_conv = e(ci_r_cl) in `row'
	qui replace lci_robust = e(ci_l_rb) in `row'
	qui replace uci_robust = e(ci_r_rb) in `row'
	qui replace bw_lo = e(h_l) in `row'
	qui replace bw_hi = e(h_r) in `row'
	qui replace nobs_orig = e(N) in `row'
	qui replace nobs_left = e(N_h_l) in `row'
	qui replace nobs_right = e(N_h_r) in `row'
	qui unique stdt if temp_in_reg==1
	qui replace ndist = r(unique) in `row'
	qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
	qui replace ymean = r(mean) in `row'
	qui replace ftag = "`ftag'" in `row'
	qui replace fals = "`fals'" in `row'

	// Residualize `yvar' for regression sample
	qui reg `yvar' `controls' `fe' if temp_in_reg==1
	cap drop y_resid_sample
	predict y_resid_sample, residuals
	
	// Residualize `yvar' for [100,500] bandwidth
	qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
	cap drop y_resid_full
	predict y_resid_full, residuals
	
	// RD plots for regression sample (variable bandwidth) 
	foreach bins in 10 20 40 {
		if `h_l'<=75 {
			local xlab = "250 275 300 325 350"
		}
		else if inrange(`h_l',76,125) {
			local xlab = "200 250 300 350 400"
		}
		else if inrange(`h_l',126,175) {
			local xlab = "150 200 250 300 350 400 450"
		}
		else {
			local xlab = ""
		}
		rdplot y_resid_sample tot_p if temp_in_reg==1, ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`y' brightness residuals", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(`xlab', labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
	}

	// RD plots for regression sample (constant bandwidth) 
	foreach bins in 10 20 40 {
		rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`y' brightness residuals", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(, labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
	}			
}

	// Store results in temp file
keep fs_step-ftag fals
dropmiss, obs force
tempfile fals2
save `fals2'



	// 2.3 11th Plan, multi habitation
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in FALSIFICATION sample: 11th Plan, multi habitation
gen in_fs_sampleFALS = vplan4==11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==0 & pop_non_zero==1 
keep if in_fs_sampleFALS==1
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==6 // Haryana has no villages in bandwidth
drop if st_code==14 // Manipur has no villages in regression
drop if st_code==32 // Kerala has no villages in regression
tab st_code, gen(STFE)	
drop STFE1 // to avoid collinearity

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)

  // Create variables to store regresson results
gen fs_step = .
gen yvar = ""
gen ifs = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen ndist = . 
gen ymean = .
gen ftag = ""
gen fals = ""
local folder = "RDROBUST plots fs lights"
local h_wide = 200

local fs_step = 13
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "mserd"
local vce = ""
local poly = 1
local fals = "11th Plan, multi hab"
local row = 0

foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

	// Define outcome variable
	local yvar = "lights_max`y'_hat"
	local title = "`y' Brightness, 11th Plan, Multiple Habitations"
	local ftag = "FALS_11th_multiH"

	// Remove controls used to project outcome variable
	local controls = "`CONTROLS'"
	if regexm("`yvar'","2007")==1 {
		local controls = subinstr("`controls'"," lights_max2005","",1)
	}
	else if regexm("`yvar'","2006")==1 {
		local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2005")==1 {
		local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2004")==1 {
		local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2003")==1 {
		local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2002")==1 {
		local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}

	// Run first-stage regresssion
	rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') bwselect(`bwmethod') p(`poly') all `vce'

	// Generate in-sample indicator and store bandwidths
	cap drop temp_in_reg
	qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
	qui count if temp_in_reg
	assert r(N)==(e(N_h_l)+e(N_h_r))
	local h_l = e(h_l)
	local h_r = e(h_r)

	// Store results
	local row = `row' + 1
	qui replace fs_step = `fs_step' in `row'
	qui replace yvar = "`yvar'" in `row'
	qui replace ifs = "`ifs'" in `row'
	qui replace control = "`controls'" in `row'
	qui replace fe = "`fe'" in `row'
	qui replace kernel = e(kernel) in `row'
	qui replace bwmethod = e(bwselect) in `row'
	qui replace vce = "`vce'" in `row'
	qui replace polynomial_order = e(p) in `row'
	qui replace beta_conv = e(tau_cl) in `row'
	qui replace beta_robust = e(tau_bc) in `row'
	qui replace se_conv = e(se_tau_cl) in `row'
	qui replace se_robust = e(se_tau_rb) in `row'
	qui replace pval_conv = e(pv_cl) in `row'
	qui replace pval_robust = e(pv_rb) in `row'
	qui replace lci_conv = e(ci_l_cl) in `row'
	qui replace uci_conv = e(ci_r_cl) in `row'
	qui replace lci_robust = e(ci_l_rb) in `row'
	qui replace uci_robust = e(ci_r_rb) in `row'
	qui replace bw_lo = e(h_l) in `row'
	qui replace bw_hi = e(h_r) in `row'
	qui replace nobs_orig = e(N) in `row'
	qui replace nobs_left = e(N_h_l) in `row'
	qui replace nobs_right = e(N_h_r) in `row'
	qui unique stdt if temp_in_reg==1
	qui replace ndist = r(unique) in `row'
	qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
	qui replace ymean = r(mean) in `row'
	qui replace ftag = "`ftag'" in `row'
	qui replace fals = "`fals'" in `row'

	// Residualize `yvar' for regression sample
	qui reg `yvar' `controls' `fe' if temp_in_reg==1
	cap drop y_resid_sample
	predict y_resid_sample, residuals
	
	// Residualize `yvar' for [100,500] bandwidth
	qui reg `yvar' `controls' `fe' if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide')
	cap drop y_resid_full
	predict y_resid_full, residuals
	
	// RD plots for regression sample (variable bandwidth) 
	foreach bins in 10 20 40 {
		if `h_l'<=75 {
			local xlab = "250 275 300 325 350"
		}
		else if inrange(`h_l',76,125) {
			local xlab = "200 250 300 350 400"
		}
		else if inrange(`h_l',126,175) {
			local xlab = "150 200 250 300 350 400 450"
		}
		else {
			local xlab = ""
		}
		rdplot y_resid_sample tot_p if temp_in_reg==1, ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`y' brightness residuals", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(`xlab', labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/lights_fs_`ftag'_`y'_inreg_`bins'.pdf", replace
	}

	// RD plots for regression sample (constant bandwidth) 
	foreach bins in 10 20 40 {
		rdplot y_resid_full tot_p if `ifs' & inrange(tot_p,299.5-`h_wide',299.5+`h_wide'), ///
			c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_wide' `h_wide') ///
			graph_options( ///
			title("`title'", color(black) size(large)) ///
			ytitle("`y' brightness residuals", size(medlarge)) ///
			xtitle("2001 village population", size(medlarge)) ///
			ylabel(,nogrid angle(0) labsize(medlarge)) ///
			xlabel(, labsize(medlarge)) ///
			graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
			legend(off))
		graph export "$results/`folder'/lights_fs_`ftag'_`y'_wide_`bins'.pdf", replace
	}			
}

	// Store results in temp file
keep fs_step-ftag fals
dropmiss, obs force
tempfile fals3
save `fals3'


	// Append falsification results to lights first stage results file
use "$results/RDROBUST_fs_lights.dta", replace
cap drop if fals!=""
append using `fals1' `fals2' `fals3'	
duplicates drop
compress
save "$results/RDROBUST_fs_lights.dta", replace

}

****************************************************************** 
****************************************************************** 

** 3. Adhoc bandwidth sensitivity analysis (manual bandwidths)
{
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in RD sample
gen in_fs_sample = vplan4<11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==1 & pop_non_zero==1 
keep if in_fs_sample==1
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==32 // Kerala, only 3 villages
tab st_code, gen(STFE)	
drop STFE1 // to avoid collinearity

	// Create district FEs (since RD robust doesn't let you pass them through)
gen stdtFE = stdt
tab stdtFE state     if inlist(stdtFE,62,63,64,70,91,180,181,202,308,309,320,490,493,494,495,504,505,508)	
replace stdtFE = 999 if inlist(stdtFE,62,63,64,70,91,180,181,202,308,309,320,490,493,494,495,504,505,508)	
	// one catch-all district FE for districts with so few in-sample villages that they break rdrobust
tab stdtFE, gen(DTFE)	
drop DTFE1 // to avoid collinearity

	// Create block groups (for clustering)
egen stdtbk = group(stdt bk_code)

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)


  // Create variables to store regresson results
gen fs_step = .
gen yvar = ""
gen ifs = ""
gen control = ""
gen fe = ""
gen kernel = ""
gen bwmethod = ""
gen vce = ""
gen polynomial_order = .
gen beta_conv = .
gen beta_robust = .
gen se_conv = .
gen se_robust = .
gen pval_conv = .
gen pval_robust = .
gen lci_conv = .
gen uci_conv = .
gen lci_robust = .
gen uci_robust = .
gen bw_lo = .
gen bw_hi = .
gen nobs_orig = .
gen nobs_left = .
gen nobs_right = .
gen ndist = . 
gen ymean = .
gen ftag = ""
local folder = "RDROBUST plots fs lights"

	// 3.1 Hard-coded bandwidths
local fs_step = 14
local ifs = "pop_mismatch20==0 & lights_diff<20"
local CONTROLS = "lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005" 
local fe = "STFE*"
local kernel = "tri"
local bwmethod = "manual"
local vce = ""
local poly = 1
local row = 0

foreach y in 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 {

	// Define outcome variable
	local yvar = "lights_max`y'_hat"
	local title = "Preferred specification: `y' Brightness"
	local ftag = "Manual bandwidth"

	// Remove controls used to project outcome variable
	local controls = "`CONTROLS'"
	if regexm("`yvar'","2007")==1 {
		local controls = subinstr("`controls'"," lights_max2005","",1)
	}
	else if regexm("`yvar'","2006")==1 {
		local controls = subinstr("`controls'"," lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2005")==1 {
		local controls = subinstr("`controls'"," lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2004")==1 {
		local controls = subinstr("`controls'"," lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2003")==1 {
		local controls = subinstr("`controls'"," lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	else if regexm("`yvar'","2002")==1 {
		local controls = subinstr("`controls'"," lights_max2000 lights_max2001 lights_max2002 lights_max2003 lights_max2004 lights_max2005","",1)
	}
	
	foreach h in 50 60 70 80 90 100 110 120 130 140 150 {
		
		// Run first-stage regresssion
		rdrobust `yvar' tot_p if `ifs', c(299.5) covs(`controls' `fe') kernel(`kernel') h(`h' `h') p(`poly') all `vce'

		// Generate in-sample indicator and store bandwidths
		cap drop temp_in_reg
		qui gen temp_in_reg = `ifs' & inrange(tot_p,299.5-e(h_l),299.5+e(h_r))
		qui count if temp_in_reg
		assert r(N)==(e(N_h_l)+e(N_h_r))
		local h_l = e(h_l)
		local h_r = e(h_r)

		// Store results
		local row = `row' + 1
		qui replace fs_step = `fs_step' in `row'
		qui replace yvar = "`yvar'" in `row'
		qui replace ifs = "`ifs'" in `row'
		qui replace control = "`controls'" in `row'
		qui replace fe = "`fe'" in `row'
		qui replace kernel = e(kernel) in `row'
		qui replace bwmethod = "manual" in `row'
		qui replace vce = "`vce'" in `row'
		qui replace polynomial_order = e(p) in `row'
		qui replace beta_conv = e(tau_cl) in `row'
		qui replace beta_robust = e(tau_bc) in `row'
		qui replace se_conv = e(se_tau_cl) in `row'
		qui replace se_robust = e(se_tau_rb) in `row'
		qui replace pval_conv = e(pv_cl) in `row'
		qui replace pval_robust = e(pv_rb) in `row'
		qui replace lci_conv = e(ci_l_cl) in `row'
		qui replace uci_conv = e(ci_r_cl) in `row'
		qui replace lci_robust = e(ci_l_rb) in `row'
		qui replace uci_robust = e(ci_r_rb) in `row'
		qui replace bw_lo = e(h_l) in `row'
		qui replace bw_hi = e(h_r) in `row'
		qui replace nobs_orig = e(N) in `row'
		qui replace nobs_left = e(N_h_l) in `row'
		qui replace nobs_right = e(N_h_r) in `row'
		qui unique stdt if temp_in_reg==1
		qui replace ndist = r(unique) in `row'
		qui sum `yvar' if temp_in_reg==1 & tot_p<299.5
		qui replace ymean = r(mean) in `row'
		qui replace ftag = "`ftag'" in `row'
	
		// Residualize `yvar' for regression sample
		qui reg `yvar' `controls' `fe' if temp_in_reg==1
		cap drop y_resid_sample
		predict y_resid_sample, residuals
		
		// RD plots for regression sample (variable bandwidth) 
		foreach bins in 10 20 40 {
			if `h_l'<=75 {
				local xlab = "250 275 300 325 350"
			}
			else if inrange(`h_l',76,125) {
				local xlab = "200 250 300 350 400"
			}
			else if inrange(`h_l',126,175) {
				local xlab = "150 200 250 300 350 400 450"
			}
			else {
				local xlab = ""
			}
			rdplot y_resid_sample tot_p if temp_in_reg==1, ///
				c(299.5) p(`poly') nbins(`bins' `bins') kernel(`kernel') h(`h_l' `h_r') ///
				graph_options( ///
				title("`title'", color(black) size(large)) ///
				ytitle("`y' brightness residuals", size(medlarge)) ///
				xtitle("2001 village population", size(medlarge)) ///
				ylabel(,nogrid angle(0) labsize(medlarge)) ///
				xlabel(`xlab', labsize(medlarge)) ///
				graphregion(color(white)) plotregion(fcolor(white)) graphregion(lcolor(white)) ///
				legend(off))
			graph export "$results/`folder'/lights_fs_manual`h'_`y'_inreg_`bins'.pdf", replace
		}
	}			
}


	// Append falsification results to lights first stage results file
keep fs_step-ftag
dropmiss, obs force
tempfile bw_sens
save `bw_sens'
use "$results/RDROBUST_fs_lights.dta", replace
append using `bw_sens'
duplicates drop
compress
save "$results/RDROBUST_fs_lights.dta", replace

}

****************************************************************** 
****************************************************************** 

** 4. Placebo test for main lights RD (2011)
{
	// Export data 
use "$panel/panel_dataset_full.dta", clear

	// Keep villages in RD sample
gen in_fs_sample = vplan4<11 & corr_state==1 & sample==1 & sample_h==1 & sing_h==1 & pop_non_zero==1 
keep if in_fs_sample==1
keep if pop_mismatch20==0 
	
	// Create state FEs (since RD robust doesn't let you pass them through)
drop if st_code==32 // Kerala, only 3 villages
tab st_code, gen(STFE)	
drop STFE1 // to avoid collinearity

	// Create lights-difference variable, to identify crazy outliers
gen lights_diff = abs(lights_max2011_hat - lights_max2001_hat)
keep if lights_diff<20

	// Keep essential variables only
keep lights_max2011_hat tot_p lights_max1998 lights_max1999 lights_max2000 lights_max2001 lights_max2002 ///
	lights_max2003 lights_max2004 lights_max2005 STFE* 

	// Save dataset for placebo tests
save "$sens/RDROBUST_placebo_test_lights.dta", replace
	
	// Run "PLACEBO_test_rdrobust_lights.do"
	
}

****************************************************************** 
****************************************************************** 







